NASA/TM-200 1 -2 1 039 1 


\ 



Nonstationary Dynamics Data Analysis With 
Wavelet-SVD Filtering 


Marty Brenner 

NASA Dryden Flight Research Center 
Edwards, California 

Dale Groutage 

Naval Surface Warfare Center (NSWC) 
Bremerton, Washington 


April 2001 



The NASA STI Program Office... in Profile 


Since its founding, NASA has been dedicated 
to the advancement of aeronautics and space 
science. The NASA Scientific and Technical 
Information (STI) Program Office plays a key 
part in helping NASA maintain this 
important role. 

The NASA STI Program Office is operated by 
Langley Research Center, the lead center for 
NASA’s scientific and technical information. 

The NASA STI Program Office provides access 
to the NASA STI Database, the largest collection 
of aeronautical and space science STI in the 
world. The Program Office is also NASA's 
institutional mechanism for disseminating the 
results of its research and development activities. 
These results are published by NASA in the 
NASA STI Report Series, which includes the 
following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant 
phase of research that present the results of 
NASA programs and include extensive data 
or theoretical analysis. Includes compilations 
of significant scientific and technical data 
and information deemed to be of continuing 
reference value. NASA’s counterpart of 
peer-reviewed formal professional papers but 
has less stringent limitations on manuscript 
length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or 
of specialized interest, e.g., quick release 
reports, working papers, and bibliographies 
that contain minimal annotation. Does not 
contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. 

Collected papers from scientific and 
technical conferences, symposia, seminars, 
or other meetings sponsored or cosponsored 
by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and mission, 
often concerned with subjects having 
substantial public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific 
and technical material pertinent to 
NASA’s mission. 

Specialized services that complement the STI 
Program Office’s diverse offerings include 
creating custom thesauri, building customized 
databases, organizing and publishing research 
results. . .even providing videos. 

For more information about the NASA STI 
Program Office, see the following: 

• Access the NASA STI Program Home Page 
at http://www.sti.nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA Access Help 
Desk at (301) 621-0134 

• Telephone the NASA Access Help Desk at 
(301)621-0390 

• Write to: 

NASA Access Help Desk 
NASA Center for AeroSpace Information 
7121 Standard Drive 
Hanover, MD 21076-1320 


t 



NAS A/TM-200 1 -2 1 039 1 



Nonstationary Dynamics Data Analysis With 
Wavelet-SVD Filtering 


Marty Brenner 

NASA Dryden Flight Research Center 
Edwards, California 

Dale Groutage 

Naval Surface Warfare Center (NSWC) 
Bremerton, Washington 


National Aeronautics and 
Space Administration 

Dryden Flight Research Center 
Edwards, California 93523-0273 


April 2001 



Available from the following: 


NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301)621-0390 


National Technical Information Service (NTIS) 

5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 487-4650 



NONSTATIONARY DYNAMICS DATA ANALYSIS 
WITH WAVELET-SVD FILTERING 


Marty Brenner* 
Aerostructures Branch 
NASA Dryden Flight Research Center 
Edwards, CA 93523-0273 


Abstract 

Nonstationary time-frequency analysis is used for identifica- 
tion and classification of aeroelastic and aeroservoelastic dy- 
namics. Time-frequency multiscale wavelet processing gen- 
erates discrete energy density distributions. The distribu- 
tions are processed using the singular value decomposition 
(SVD). Discrete density functions derived from the SVD gen- 
erate moments that detect the principal features in the data. 
The SVD standard basis vectors are applied and then com- 
pared with a transformed- SVD , or TSVD, which reduces the 
number of features into more compact energy density concen- 
trations. Finally, from the feature extraction, wavelet-based 
modal parameter estimation is applied. 

The primary objective is the automation of time-frequency 
analysis with modal system identification. The contribution 
is a more general approach in which distinct analysis tools are 
merged into a unified procedure for linear and nonlinear data 
analysis. This method is first applied to aeroelastic pitch- 
plunge wing section models. Instability is detected in the 
linear system, and nonlinear dynamics are observed from the 
time-frequency map and parameter estimates of the nonlinear 
system. Aeroelastic and aeroservoelastic flight data from the 
DAST (drone for aerodynamic and structural testing) and 
F18 aircraft are also investigated and comparisons made be- 
tween the SVD and TSVD results. Input-output data is used 
to show that this process is an efficient and reliable tool for 
automated on-line analysis. 


Nomenclature 

a wavelet scale 

DAST drone for aerodynamic and 

structural testing 

DWT discrete wavelet transform 
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Fi 

time-frequency feature descriptor 

9 

wavelet basis function 

j 

imaginary unit 

i, k, l 

integer indices 

m, n 

number of time and frequency points 

p,q 

moment orders 

r,s,t 

integer indices 

R 

rank of time-frequency map 

SRA 

Systems Research Aircraft 

SVD 

singular value decomposition 

TSVD 

transformed- SVD 

w. 

continuous wavelet transform 

T 

wavelet translation time 

UJ 

radian frequency 


wavelet peak frequency 

* 

complex conjugate transpose 

<> 

result of moment calculation 


Introduction 

Many aerospace disciplines face the task of analyzing 
nonstationary signals in which the frequency content 
changes with time. Most physical phenomena from 
acoustics, aerodynamics, thermodynamics, structural 
dynamics, propulsion, and controls are naturally time- 
varying events with frequency variations, transients, and 
complex harmonic interactions. Requirements for adap- 
tive nonlinear procedures are now becoming more essen- 
tial components of modern data analysis of aerospace 
systems in concordance with numerous other engineer- 
ing and scientific applications. 

A fundamental objective in data analysis of physical sys- 
tems is to obtain accurate representation of the dynamics 
for a particular analysis. Optimal and relevant represen- 
tations are desired for efficiency and consistency between 
the data decomposition and the physical system I 4, 28, 31 ^. 
Methods for achieving this goal are often based on en- 
tropy or other statistical measures for best-basis decom- 
positions and optimal performance ^ 22, 42 ^. 
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Parsimonious representation of data has many such im- 
plications in denoising, compression, information re- 
trieval, detection, identification, classification, and pat- 
tern recognition This paper emphasizes applications 
in dynamics analysis of aeroelastic and aeroservoelastic 
systems from test data. Dominant and important fea- 
tures in the data are extracted to provide automated 
information retrieval for system identification. The data 
analysis becomes an integral part of the system identifi- 
cation procedure for tracking modal stability estimates 
in a nonstationary setting. 

Time-frequency analysis provides a powerful tool for the 
analysis of nonstationary signals t 18 l Applications have 
been demonstrated in music signal analysis machin- 
ery diagnosis [26 ', damage detection 1401 , seismic moni- 
toring t 19 l, medical signal analysis t 41 >, and acoustic and 
speech processing t 32 L Transfer functions and modal 
parameter estimates derived from time-frequency rep- 
resentations have been applied to estimate state-space 
aeroservoelastic models l 2, 3| 10 ). 

Time- frequency structure is revealed by quantifying the 
distribution of signal energy as a joint function of time 
and frequency. Localization of the energy density de- 
scribes energy density concentrations at specific loca- 
tions in the time-frequency plane. Dominant and im- 
portant concentrations need to be accurately discrimi- 
nated to provide descriptors relating to the location in 
time, time duration, frequency location, and local band- 
width of principal energy density l n » 17> 431 . Minimizing 
the number of these descriptors while preserving salient 
information from the energy distribution for modal esti- 
mation is desired for efficient dynamics analysis. 

Many approaches have been proposed to extract a sig- 
nal from noise in the time-frequency representation I 33 l 
One approach performs a masking operation on the 
time-frequency map, then separates out desirable fea- 
tures Often the regions of the signal are well- 

concentrated relative to the widely-distributed noise so 
this method works well in general. However, it requires 
user intervention unless regions are well-defined apriori. 
For this reason, an automated filtering procedure is pro- 
posed using the singular value decomposition (SVD) of 
the time-frequency data I 29, 43 K 

The singular value decomposition is ubiquitous in the sig- 
nal processing, identification, and control fields as a tool 
for rank-revealing algorithms, model reduction, subspace 
detection, and information retrieval l 1 * 20) . The SVD is 
applied in this paper to derive density functions for gen- 
erating moments i 14] from the joint energy density time- 
frequency distribution. Moments relate to the principal 
signatures of the nonstationary signals ^ 27 1. In the con- 
text of subspace estimation and information retrieval it 
is a data reduction for model identification. 


The SVD representation tends to be sparse relative to 
an entire time- frequency map. The orthonormality of 
the SVD singular vectors for density extraction allows 
independent decomposition of time and frequency con- 
tent in the data, thereby avoiding joint moment consid- 
erations. A new transformation of the SVD basis vec- 
tors, called the transformed-SVD (TSVD) ^ 15 1, is also 
used to produce more concentrated descriptions of en- 
ergy density for multicomponent signals containing sim- 
ilar components distributed in the time-frequency plane. 
For separation of components requiring very fine time or 
frequency resolution, the TSVD may be advantageous. 

Multiresolution wavelet signal processing has shown 
promise for studying time-frequency characteristics of 
signals by decomposing data into cells with properties 
of scale and frequency concentrated in time. These cells 
form a type of tiling and consist of Gaussian-windowed 
sinusoid basis functions (atoms), also known as Morlet 
wavelets, creating an atomic multiscale decomposition 
in a filter bank structure I 44] . Competing requirements 
of time and frequency resolution, subject to the uncer- 
tainty principle I 28 *, is accomplished with a combination 
of dyadic multiscale decomposition, compact orthogonal- 
ity, and harmonic wavelet properties f 23, 24 l 

This paper augments time-frequency multiscale wavelet 
processing with SVD filtering and wavelet-based modal 
parameter estimation. The contribution is a more gen- 
eral approach in which distinct analysis tools are merged 
into a unified procedure: 

• multiresolution analysis with wavelet decomposi- 
tion at multiple scales 

• SVD techniques for information retrieval 

• feature extraction from moments of densities 

• modal parameter estimation from the complex 
Morlet basis functions 

This automated nonlinear filtering procedure is first 
applied to linear and nonlinear two degree-of-freedom 
pitch-plunge wing section testbed models. Nonlinear dy- 
namics are detected from the time-frequency map of the 
nonlinear system. Aeroelastic dynamics from the experi- 
mental flutter flight test drone, DAST (drone for aerody- 
namic and structural testing) , is used to show how the in- 
put time-frequency signature can be used for automating 
the input-output data analysis. Data from a flutter en- 
counter shows the capability for detecting an approach- 
ing instability from actual flight test data. Aeroservoe- 
lastic flight data from the F18 Systems Research Air- 
craft (SRA) is also investigated and comparisons made 
between the SVD and TSVD results. 
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Wavelet-SVD Method 

This section describes the method to construct time- 
frequency signed representations from a Morlet wavelet 
decomposition, then to use this energy distribution map 
to extract principal features from the SVD. Density func- 
tions derived from the SVD singular vectors generate 
spectral and temporal moments relating to features of 
the original nonstationary time series process. These 
moments are used for subsequent wavelet-based modal 
parameter estimation and system identification. 


Multiscale transforms are used to exploit multiresolution 
analysis, or redundant representations of a signal on mul- 
tiple frequency bands. Nonorthonormal Morlet wavelets 
are approximated with harmonic-like discretizations on 
wavelet scales corresponding to subharmonics represent- 
ing overlapping frequency bands per octave (where each 
octave is a doubling of resolution, or a scale twice as 
small). These wavelets form a nonorthogonal redundant 
basis for the signal space ^ . Adjustment to satisfy the 
competing requirements of time and frequency resolution 
is accomplished with a combination °^°™P act 
onal and harmonic wavelet properties I 23, 24 . 


Time -Frequency Wavelet Decomposition 

Much progress has been made in the development of 
time-frequency distributions ^ . For a time- frequency 
distribution Q to be interpreted as a joint time-frequency 
energy density, it must be nonnegative and satisfy the 
correct time and frequency marginals for all time and 
frequency 


where X(f) = fZ exp(-j2ir/t)* is the Fourier 
transform of the signal x(t). Marginals \x(t)\ 2 and 
\X(f)\ 2 are the energy densities of time and frequency 
which are commonly called the instantaneous power and 
the energy density spectrum, respectively. Algorithms 
for constructing proper positive distributions satisfying 
these properties have recently been developed [16 ’ 37 . 
The joint energy density specifies concentrations of en- 
ergy in time and frequency. 

Feature extraction from joint distributions proceeds from 
the joint time- frequency moments of signal x(t) (assumed 
to have unit energy) given by 

<it> 

p,q = 1,2,3, ••• (1) 

which defines the temporal and spectral moments 

/ OO 

t p \x(t)\ 2 dt 

-OO 

/ OO 

f q \x(f)\ 2 df. ( 2 ) 

-OO 


/ OO 

Q(tJ)df = \x(t)\ 2 

-OO 


In the present application, the time-frequency distribu- 
tion, Q{t,f), is constructed from a multiscale wavelet 
transform l 3 l Therefore, the distribution is not the strict 
form of a joint energy density since the requirement on 
marginals is not enforced. The wavelet transform tends 
to spread the energy over time and frequency, yet the 
wavelet decomposition is essential for identification of 
modal dynamics as will be shown. This energy density 
approximation is therefore relevant for this problem. 


A discrete wavelet transform (DWT) is derived from 
the wavelet basis to get a multiresolution analysis of 
the sampled continuous Morlet transform. Fast algo- 
rithms are realized with a dyadic multiscale decomposi- 
tion. Therefore, the DWT is implemented as a dyadic 
filter bank covering a pre-defined range of frequencies 
with corresponding number of frequency bands per oc- 
tave (voices/octave) I 44 '. Voices can be viewed as frac- 
tional dilations of a single wavelet at a particular scale. 
Approximate orthogonality is imposed in the multires- 
olution representation of the Morlet filters. Multiscale 
Morlet DWTs provide efficient and flexible means for 
analysis of nonstationary data with adjustable frequency 
resolution versus time localization. 


The continuous wavelet transform of signal x(t) over the 
time-scale (r, a) plane is represented as 


Wg(T,a) = ^J *Ws* ( Sr ) dt 

where scale parameter a is proportional to the duration 
and inversely proportional to the peak frequency w 0 of 
the complex Morlet wavelet 


g{t) = 


s/2n 




e 2 e 


The spectrum of a dilated and translated Morlet wavelet 

G„jLj)=e- {au - u ° ) 2 e juT 


reaches a maximum value at a = “f. The continuous 
Morlet basis functions are discretized for a DWT filter 
bank decomposition. Signal power is preserved by en- 
forcing the following identity with the DWT {C w is an 
admissibility constant) 

from which the instantaneous power is expressed in terms 
of the wavelet time variable, r, as 

= XT |W f (r,.)|»£. 

JO 

This time-scale decomposition of data is often called a 
scalogram, or the energy density |W 3 (r,a)| 2 of the signal 
over the (r, a) plane I 28 ' . 
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As an example of a bank of discretized finite impulse 
response Morlet filters, figure 1 shows an array of wavelet 
filters for a number of voices per octave. These filters are 
adjusted on the plot to line up the peak magnitudes at 
the same locations in time. Center frequency increases 
with voice number, and the filters are normalized to have 
the same magnitude in the frequency domain at their 
respective peak frequencies. 


FIR Discrete Wavelet Filters 



Filter Length, sec ® ® Voice Number 


Figure 1: Wavelet filters for a number of voices per octave. 

Morlet wavelets are being used to create the power distri- 
bution and to estimate the modal parameters. Therefore, 
an implicit filtering process is being performed indepen- 
dent of the explicit procedure of feature extraction and 
noise removal using the SVD and TSVD. The wavelet ba- 
sis representation of the signal is a projection subspace 
from which modal parameters are derived. 

SVD Filtering 

Numerous wavelet thresholding techniques exist based 
on statistical measures of significance for denoising. 
These include entropy measures from information and 
communication theory. They are generally based on reg- 
ularity, smoothness, or a low-pass character of the un- 
derlying signal relative to noise, and the compromise is 
between goodness-of-fit and smoothness. Because this 
approach is somewhat subjective by nature, a more de- 
terministic approach is preferable. In aeroelasticity and 
aeroservoelasticity, the statistical criteria are generally 
not applicable since the dynamics can be of significant 
bandwidth and include high-frequency transients. 

Assume the time-frequency distribution is a positive dis- 
tribution in the form of a scalogram, spectrogram, or 
any other representation such as the Wigner-Ville dis- 
tribution. In the present application the basis func- 
tions are the discretized complex Morlet wavelets, so 
the scalogram magnitude is appropriate. Discretized 
wavelets form the multiresolution analysis which can be 
exploited at different scales for information extraction 
by the SVD. 


The SVD spectrum encodes the features of time- 
bandwidth product, frequency-time dependence, and 
number and location of signal components. It is invari- 
ant to shifts of the signal in time or frequency [29 ’ 43) . 
While this SVD method is described for the entire time- 
frequency distribution, it could also be used to analyze 
the different scales independently or jointly. 

Consider the distribution Q as an m-by-n discrete matrix 
A such that A(kJ) = Q{tk,fi ), where k - 1,2 
and / = 1, 2, • • • , n for m time points and n discrete fre- 
quencies. This matrix can always be decomposed with 
the SVD into a set of basis matrices A x with correspond- 
ing singular values, or weights, a x 

R 

A = ^ a,A t ; At = u t v-. (3) 

i - 1 

R is the rank, the set of constant weights a x are ordered 
such that G\ > CT2 > • " > VR > and each A i ( or a i) 
corresponds to singular vectors ti; and as shown. Time 
and frequency aspects of the A t are independent since 
m and are orthonormal. The maximum value of the 
rank R can only be as large as the minimum of the row- 
column (m-n) dimension of A. The objective is to get a 
reduced set of energy concentrations defined by the A x 
to generate time-frequency features based on the relative 
magnitudes of the associated a x . Note that the rank R 
can never exceed the number of discrete frequencies n, 
i.e., R < n, since there will always be more time points 
than discrete frequencies. 

The matrix A is used with the weights to extract high- 
lights from the wavelet map by characterizing the A x 
with the joint distribution of equation 1 . For each a x , a 
corresponding region in the distribution is desired (At) 
which is weighted in value by the magnitude of a x re- 
lating to its respective information content. From the 
original joint moments of equation 1, the SVD separates 
the distribution into R primary moments. This decom- 
position correlates directly with independent time and 
spectral moments expressed in equation 2 through inde- 
pendent singular vectors tq and V( in equation 3, so joint 
moments need not be considered 

A unit-energy positive distribution A matrix does not 
guarantee that each Aj will be element-by-element pos- 
itive, which is necessary for a proper density function. 
Matrices A x are constructed from the A x by squaring 
each matrix element to create proper densities for mo- 
ment calculations. Each A* is formed from the element- 
by-element square of the associated singular vectors u x 
and Vi , and the corresponding squared vectors are de- 
noted by Ui and Vi . The new basis matrices are 

Ai — u x v* > 0. 

Density functions u x and u* compose proper density func- 
tions Ai by being positive and orthonormal, and by be- 
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ing independent, they are also the time and frequency 
marginals of the respective Moments of the den- 
sity functions constructed from the singular time and 
frequency vectors are used to formulate features. Tem- 
poral and spectral moments corresponding to equation 2 
are calculated as 

m ™ 

< >,= y>) p M*); < f q >•= 

*=i ' =i 

First moments (p = q = 1) estimate time-correlated in- 
stantaneous frequency, and the second moments (p — 
q = 2) determine the time duration and instantaneous 
bandwidth for each feature in the distribution. A feature 
descriptor F{ is therefore defined by 

Fi — ^<5*, /ij ti) 

= <t>i,< f >i, yj<t 2 >i -<?, \j< T 2 >> ~fPj ■ 

Each feature descriptor Fi for each .4, is a 5-element vec- 
tor which includes the normalized singular value_ weight 
di, time location U, instantaneous frequency /*, time 
duration U, and instantaneous bandwidth fi, of the fea- 
ture. Features are d r scaled rectangle regions in the time- 
frequency map. They are used as windows on the original 
scalogram to extract regions for further analysis. 

The number of desired features should preferably be as 
small as possible while preserving relevant dynamics. 
The separation property of the SVD helps maintain im- 
portant information with a minimal number of features. 
Generalizing the same process to higher-order moments 
gives instantaneous skew and kurtosis !6 - which may re- 
veal data anomalies, asymmetries, and nonlinearity. 


Simulated Four Sine-Burst Signal 



Time, sec 


Figure 2: Four sine-bursts of equal amplitude. 

A simple example of a simulated analysis with SVD filter- 
ing uses the signal in figure 2. There are four sine-bursts 
at 15, 10, 5, and 25Hz, respectively. A wavelet decompo- 
sition and SVD filtering procedure is performed on this 


signal to produce the plots in figure 3. The top plot rep- 
resents an energy density distribution contour (plotted 
in two dimensions) of the original wavelet coefficients. 
The middle plot represents the sum of all rectangular 
features Fi determined from the SVD filter. The bot- 
tom plot accentuates the middle plot using a log-scale 
(dB) for the out-of-plane contour dimension. A log-scale 
region will cover at least as much as the linear scale con- 
tour of the middle plot, and all contours are normalized 
to their maximum magnitude. The middle and bottom 
plots are identical except for the contour scale. 


7j me _F r equency Energy Density Distributions 



Figure 3: Contour plots from original wavelet decom 
position (top), SVD-filtered (middle), and SVD-filtered 
(bottom, log-scale) scalograms of four sine-burst signal. 

Note the more compact representation from SVD filter- 
ing (middle plot) compared to the original decomposi- 
tion. For an automated analysis, we choose the order 
of the distribution to be the full rank R. Of the to- 
tal 300,000 elements in the distribution, the rank in this 
case reduces the number of features to R = 300, resulting 
in a dramatic data reduction for model identification. 

Transformed- SVD 

A potential problem with the standard SVD filter stems 
from the information contained in the individual Ai 
terms of the decomposition. The principal information 
is ordered according to values of <Xi, but often there are 
multiple time-frequency regions with similar energy den- 
sity magnitudes which will contribute equally to each 
Ai . Moments calculated from the singular time and fre- 
quency vectors contain simultaneous contributions from 
distributed regions since the vectors span the entire 
scalogram. Information is not concentrated because the 
SVD cannot separate energy density concentration in the 
individual Ai so the features consist of linear combina- 
tions of many singular vectors. 
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Feature extraction would benefit from more localized sig- 
nal content in each A{ for two reasons. First, for a fixed 
number of features, fewer terms in the SVD would en- 
hance specific localized content. Second, for a fixed num- 
ber of terms, higher fidelity features would be produced. 
Therefore, what is desired is a rotation of the singular 
basis to minimize the number of basis vectors to inter- 
pret each density concentration. Ideally, each Ai term 
should relate to a single localized concentration. 

This problem cam be posed as follows: find an orthonor- 
mal transformation for the U{ and v* vectors such that 
the means are maximized in the new basis. The attempt 
is to concentrate the densities into smallest regions in the 
plane. For unknown transformation C, and orthonormal- 
ity condition CC * = C*C = /, the means are expressed 
in terms of the Ci* coefficients. Orthonormality implies 
that Vi( r ) = 1> so the means of the transformed 

vectors y* = Ylk c '^ Uk are the first moments < y r > 
given by 

<Vr >i= EHi ryf(r) 

= r £?=i c., iU . (r) £f =1 c t<i u t (r) 

= cJMc i 

In this quadratic form, the unique solution for maximiz- 
ing the means of the yi vectors is achieved when the c i are 
the eigenvectors of the M matrix. Similarly, am orthonor- 
mal matrix D is found independently for the vectors Vi by 
maximizing the means of vectors d t} kVk • Com- 

pared to the standard SVD, the transformed-SVD is as 
follows 

A = USV *; U = YC •; V = XD* 

and the proper substitutions made to get the TSVD 

A = ( YC')S{XD*Y = Y{C*SD)X* = YZX\ 

Now that instead of just a sequence of singular values a* 
from the S matrix of the SVD, there exists a full singular 
matrix Z = C*SD of terms to weight the new singular 
vectors X and Y in the TSVD. These terms have the 
same function as the terms in S, but instead of corre- 
spondence to linear combinations of concentrations in 
Ui and Vi , the Z-elements, z**, correlate to transformed 
singular vectors yi and Obviously there are signif- 
icantly more (m-by-n) of the than the Oi ( R < n). 
The singular weights in Z are ordered and the new TSVD 
decomposition for the first R of the ik- components takes 
the form (where now B ^ A, compare to equation 3) 

R 

B — ^ ^ ZikBi] Bi — y%x k * (4) 

1=1 

Returning to the signal from figure 2, a transformed- 
SVD analysis is performed and plotted in figure 4 for 


Time-Frequency Energy Density Distributions 
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Figure 4: Contour plots from original wavelet decom- 
position (top), TSVD-filtered (middle), and TSVD-filtered 
(bottom, log-scale) scalograms of four sine-burst signal. 

comparison with the SVD result in figure 3. The plots 
in figure 4 correspond to those in figure 3, so the top 
plot is identical. The middle TSVD contour plot in fig- 
ure 4 shows an enhanced version of the plot in figure 3 
by extension of the time-range estimate of three of the 
sine-bursts. Note the longer feature window based on 
the computed Fi compared to the original distribution 
contour. The burst at 5Hz (7 sec) is not distinguishable 
in the TSVD middle plot, but in the remarkably similar 
bottom log-scale plot it reappears. This burst at 5Hz 
has the least total energy content. 

Observe that in the bottom two plots, the TSVD has re- 
markably more detail with the same number (R = 300) 
of singular vectors. Information extraction is devoted to 
the regions of the sine-bursts, and not spread over unnec- 
essary areas of the time-frequency plane as in the case of 
the SVD in the bottom plot of figure 3. The more com- 
pact representation from the TSVD will improve further 
analysis by concentrating on the important regions of 
interest for modal identification. 

It was already shown in a previous study ^ 15 1 how the 
TSVD separates and improves the quality of the time- 
frequency features with a small set of descriptors for clas- 
sification purposes. The present application for flight 
data analysis emphasizes the filtering aspects of the 
TSVD versus the SVD. The purpose is not to find a min- 
imum number of descriptors, but to get the highest qual- 
ity filtered output for a fixed number of features deter- 
mined by the rank, R, of the discretized time- frequency 
map. This will allow an automated procedure for pre- 
processing with a time-frequency filter to augment modal 
parameter estimation and system identification. 
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SVD approaches are attractive from a systems analysis 
perspective. Prevalence in model reduction, subspace 
identification, and optimal filtering point to its signif- 
icance as a tool for information retrieval. Some of its 
power lies in the geometry of the singular vectors as a 
distance measure between subspaces l 7 !. Also, for on-line 
implementations, popular SVD recursive procedures al- 
low for adaptive signal processing by using fast updating- 
downdating schemes t 1 - 9 '. The idea is to perform the 
wavelet decompositions, update the time-varying distri- 
butions, and apply on-line procedures for the analysis. 

Paramptric Modal Estimation 

Modal parameters can be estimated with wavelets by 
analysis of the system impulse response I 3 ' 36) (see ap- 
pendix). The DWT of a signal using the complex Morlet 
wavelet is a complex-valued matrix whose modulus and 
phase are related to impulse response parameters. Linear 
phase variation of the DWT estimates the instantaneous 
frequency. Wavelet modulus decay is used similarly to 
derive decay rate for corresponding modal damping es- 
timates assuming the instantaneous frequencies corre- 
spond to modal dynamics. In the current application, 
this procedure is applied at every time point assuming 
at each instant that the response is a sum of multiple 
degree-of-freedom impulse responses. 



Time, sec Time, sec 

Figure 5: Instantaneous frequency estimates (left) and 
detected dominant spectrum (right) using the SVD. 

Instantaneous frequency estimation of the four sine-burst 
signal of figure 2 is illustrated in figures 5 and 6 using 
the SVD and TSVD, respectively. The plots show fre- 
quency estimates over the corresponding detected period 
of time from the DWT. Right plots titled Dominant Fre- 
quencies in each figure are subsets of the those titled Nat- 
ural Frequency. Left plots designate raw estimates from 
the respective filtered regions, and right plots refer to the 
spectrum detected by using a threshold parameter on the 
wavelet coefficients. This threshold is 0.5 of the normal- 
ized absolute magnitude of the scalogram. Threshold 
estimates are referred to in this report as dominant. 


In the SVD analysis of figure 5, numerous frequencies are 
detected from the extracted regions of figure 3, but of 
those only two dominant frequencies at 15Hz (3 sec) and 
25Hz (9 sec) in the right plot are estimated based on this 
threshold. The congested left plot indicates that many 
detected frequencies are redundant or spurious estimates 
from large relatively insignificant regions of the time- 
frequency map (bottom plot of figure 3). Conversely, 
in figure 6 for the TSVD, the left and right plots are 
identical and all the frequencies are detected as being 
dominant. This demonstrates that the TSVD is detect- 
ing only that part of the spectrum that is significant in 
the data without any threshold criteria. 
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Figure 6: Instantaneous frequency estimates (left) and 
detected dominant spectrum (right) using the TSVD. 


Aeroelastic Tes tbed Model 

Previous studies have investigated applying wavelets 
for analysis of structural and aeroelastic nonlinear sys- 
tems t 25 ' 30 ’ 451 . SVDs have been used to extract the min- 
imum embedding dimension from noisy nonlinear sys- 
tems 1 39 1. In this section the wavelet-SVD approaches 
are applied to investigate linear and nonlinear reponses 
for a simple aeroelastic system. An aeroelastic testbed 
has been developed at Texas A&M University for flut- 
ter research using a prototypical aeroelastic wing sec- 
tion I 12 L This quasisteady pitch-plunge system is de- 
scribed by two-degree-of-freedom aeroelastic equations 
of motion. A control surface configuration was used for 
closed-loop control l 21 '. 

Nonlinearity is introduced to the system dynamics 
through the stiffness associated with pitch. This stiff- 
ness is described by a nonlinear polynomial function of 
the pitch angle. Such structural nonlinearities occur in 
physical aeroelastic systems and have been investigated 
to determine their effect on inducing limit cycle oscilla- 
tions. The two parameters that determine the response 
of the wing and the onset of limit cycle oscillations, are 
the elastic axis location and the freestream velocity. The 
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Linear Aeroelastic Testbed Simulation 
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Figure 7: Chirp control command input (top) and pitch 
response (bottom) from linear aeroelastic testbed model. 

stiffness functions associated with the pitch degree of 
freedom are chosen here to represent a linear spring and 
a nonlinear hardening spring. 

Wavelet Coefficients 
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Wavelet Coefficients 
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Figure 8: Modal parameter estimates from testbed 
linear simulation using SVD (top) and TSVD (bottom). 

A classical flutter instability is encountered for the lin- 
earized dynamics near a speed of 12 m/s at about 
2Hz t 25 '. Figure 7 represents a control input chirp of 


0-5Hz (top) and the unstable pitch response at 12 m/s 
(bottom). Both SVD and TSVD results are presented in 
figure 8. These plots represent the estimates and corre- 
sponding wavelet coefficient magnitudes from which the 
estimates were derived. Magnitudes of the coefficients 
determine a relative level of significance, or measure of 
observability. The upper plot from the SVD shows the 
instability detected near 2Hz. The array of zero damping 
values at frequencies below the unstable mode indicate 
an output tracking of the input (disguised by the magni- 
tude in figure 7) until an instability is detected at 1.5Hz 
near the pitch mode natural frequency. TSVD estimates 
of the same instability are more accurate near 2Hz t 25 l 
with much fewer spurious estimates. 


Nonlinear Aeroelastic Testbed Response 



Figure 9: Pitch response from nonlinear aeroelastic 
testbed model chirp command input. 
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Figure 10: Contour plots from original wavelet decom- 
position (top) and SVD-filtered (bottom, log-scale) 
scalograms of nonlinear aeroelastic testbed simulation. 


A nonlinear response from the same chirp input also at 
12 m/s with a nonlinear hardening spring in the pitch 
degree of freedom is shown in figure 9. The correspond- 
ing wavelet scalogram and SVD-filtered contour plot are 
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shown in figure 10. Note the scattering of the output 
spectrum from a well-defined chirp input, thus indicating 
the nonlinearity . There is clear indication of higher 
harmonics from the primary 2.5Hz mode at higher fre- 
quencies near 7.5Hz and 12.5Hz (the third and fifth har- 
monics). This result is also evident from the correspond- 
ing frequency and damping estimates in figure 1 1 . From 
the multitude of SVD raw estimates, this plot shows only 
the dominant coefficients using the threshold discussed 
previously. At any time during the response, at least two 
of three of the resonances are evident. 


Dominant Wavelet Coefficients 



Featured TSVD Distribution (dB) 



Figure 12: Contour plot of TSVD-filtered scalograms 
(log-scale) of nonlinear testbed simulation response. 


More precise tracking of the modes from the TSVD is 
evident from the filtered contour plot of figure 12. This 
plot corresponds to the previous contours from figure 10. 
High definition frequency estimation compared to the 
SVD is clearly seen, and results in the raw estimates 
plotted in figure 13 being identical to the TSVD domi- 
nant estimates (not shown). 


Wavelet Coefficients 



Figure 13: Modal parameter estimates from nonlinear 
aeroelastic testbed simulation using TSVD filtering. 


PAST Flutter Data Analysis 

For an automated data analysis procedure to be effec- 
tive with a known input signal, it is advantageous to use 
the time-varying input spectrum for tracking output re- 
sponses. Linear dynamics may then be discerned from 
unmodeled dynamics and nonlinearity outside of the in- 
put frequency range. For this purpose, using actual air- 
craft data, wingtip accelerometer data is analyzed from 
the NASA DAST l 13) (drone for aerodynamic and struc- 
tural testing, figure 14). Part of the DAST program in 
1980 was to pursue investigations on a drone equipped 
with a flutter suppression system to enable flight beyond 
the open-loop stability boundary. Wingtip accelerome- 
ter response data was acquired with 10-40Hz logarithmic 
chirps and 20Hz sinusoidal doublets into the aileron con- 
trol surface. 



Figure 14: NASA DAST vehicle in flight. 


From the time-frequency path of the input signal, a mask 
is created by SVD or TSVD filtering to be applied to 
the accelerometer response. An aileron chirp input (top) 
and wingtip accelerometer response (bottom) are shown 
in figure 15 at Mach 0.8 about a minute before the drone 
encountered an aeroservoelastic instability ■ 
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Wavelet Coefficients 


DAST Aileron Chirp and Wingtip Response 
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Figure 15: Aileron chirp (top) and wingtip accelerometer 
response (bottom) from DAST at Mach 0.8, 15,000 feet. 


Figure 16 shows contour plots of the aileron chirp input. 
The SVD-filtered envelope (bottom plot) bounds the in- 
put more than adequately. Also, shown in figure 17, 
the estimated modal frequency and damping show the 
primary wing bending mode near 20Hz, but also other 
minor harmonics of much smaller magnitude near 40Hz 
and 60Hz. Note that in the SVD-filtered contour plot of 
figure 16, the spread up to 60Hz occurs when the primary 
20Hz mode is excited by the chirp input (6-7 sec). 


Time- Frequency Energy Density Distributions 
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Figure 17: Modal parameter estimates from DAST 
wingtip acclerometer from chirp using SVD filtering. 


frequency window. This property of the TSVD may be 
used to separate linear output response from unmodeled 
dynamics and nonlinearity. 

The SVD filter generally encompasses a much larger 
region than the input alone, with the higher contours 
more adjacent to the dominant energy spectrum. While 
not as precise a filter as the TSVD, it detects a time- 
frequency envelope that thoroughly bounds the input. 
A TSVD filter concentrates more singular vector contri- 
butions (moments) into smaller areas than the SVD fil- 
ter. The TSVD may actually miss some input dynamics 
from its moments for a reasonable order approximation 
comparable to the SVD. A higher order approximation is 
possible with the TSVD as mentioned in the comments 
above equation 4, but this paper only addresses equal- 
order approximations between the SVD and TSVD. 


Featured TSVD Distribution (dB) 



Figure 16: Contour plots from original wavelet decom- 
position (top) and SVD-filtered (bottom, log-scale) 
scalograms of DAST aileron chirp input. 

In contrast, the TSVD-filtered contour of figure 18 dis- 
criminates discrete frequency bins directly along the 
time- varying chirp input. Dynamics outside of the time- 
frequency path of the input will not be observed. Esti- 
mates will only be derived locally along the input time- 


Figure 18: Contour plot of TSVD-filtered scalogram 
(log-scale) of DAST aileron chirp input. 

Wavelet reconstructions are created from the coefficients 
of the time-frequency map that are retained by filtering, 
masking, or thresholding. In a direct filtering or masking 
application, the quality of the wavelet reconstruction of 
a time response depends on the time-frequency window 
used on the response. Good results are obtained when 
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Reconstructed Wingtip Response 



Figure 19: Reconstructed SVD-filtered (top) and TSVD 
filtered (bottom) DAST wingtip accelerometer responses 
from DAST at Mach 0.8, 15,000 feet (see figure 15). 


all the significant dynamics are retained. For instance, 
with no filtering or masking, the reconstructed response 
is simply a projection onto the wavelet basis functions. 
This signal generally resembles the original signal be- 
cause of the good reconstruction properties from the re- 
dundancy in the multiscale wavelet transform. 


DAST Wingtip Accelerometer Response 



Figure 20: Wingtip accelerometer response from DAST 
aileron-pulsed input at Mach 0.825, 15,000 feet. 


the TSVD. Both reconstructions are modulated modal 
responses, but the TSVD consists of narrower bandpass 
filters in its modulation. 




Figure 21: Contour plots from original wavelet decom- 
position (top) and TSVD-filtered (bottom, log-scale) 
scalograms of DAST wingtip accelerometer response. 


Reconstructed Wingtip Accelerometer Response 



Figure 22: Reconstructed TSVD-filtered DAST wingtip 
accelerometer response at Mach 0.825, 15,000 feet. 


To demonstrate the estimation procedure for an actual 
aircraft instability, figure 20 represents the wingtip ac- 
celerometer response at Mach 0.825 immediately before 
the closed-loop instability. During the last three sine 
pulses the TSVD filter results are shown in the the con- 
tour plots of figure 21. The SVD filter has difficulty 
making a precise determination of the pulse energy den- 
sity map (not shown). Hacking the pulse responses is 
complete and accurate with the TSVD filter. 


For the wingtip accelerometer response of figure 15, com- 
parisons are made between the SVD (top) and TSVD- 
filtered (bottom) reconstructions in figure 19. The SVD 
filter covers more of the time-frequency plane, and this 
shows up in the reconstructed response as compared to 


Figure 22 represents a reconstructed TSVD-filtered ver- 
sion of the signal in figure 20. From the coefficients of 
the filter shown in figure 21, the reconstructed time re- 
sponse resembles an impulse response of the mode near 
20Hz. Estimated modal frequency and damping from the 
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TSVD in figure 23 shows clearly the 20Hz wing bending 
mode and reduction in modal damping with time. 

Dominant Frequencies Dominant Dampings 



Figure 23: TSVD estimated modal frequencies (left) and 
damping ratios (right) from DAST wingtip accelerometer. 


F18 Aeroservoelastic Data Analysis 


Time-Frequency Energy Density Distributions 
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Figure 25: Contour plots from original wavelet decom- 
position (top) and SVD-filtered (bottom, log-scale) 
scalograms of F18-SRA wingtip accelerometer response. 


Digital excitation signals were implemented on an the 
NASA F18 Systems Research Aircraft (SRA) for a base- 
line model update and validation effort to support the 
F18 Active Aeroelastic Wing t 35 l program. Multisine sig- 
nals were added to aileron, stabilator, and rudder actu- 
ator commands. In this example, the differential aileron 
command is used for excitation at Mach 0.9, 10, 000 feet. 
Wingtip and lateral fuselage accelerometer responses are 
used to demonstrate the SVD and TSVD filtering process 
with a discontinuous input in the time- frequency plane. 
Multisine waveforms have a chirp-like quality, but with a 
burst of high-frequency chirp at the end of the signal ^ 2 1 . 


F18-SRA Wingtip Accelerometer Response 



Featured TSVD Distribution (dB) 



Time, sec 

Figure 26: Contour plot of TSVD-filtered scalogram of 
F18-SRA multisine wingtip accelerometer response. 

SVD filter (bottom plot) attempts to bound the entire 
time-frequency response, as in the other cases, whereas 
the TSVD-filtered output in figure 26 distinguishes fine 
frequency divisions over short time intervals. The 13Hz 
and 7Hz lines near 20 sec are very distinct and precisely 
tuned to the specific frequencies. 


Figure 24: F18-SRA wingtip accelerometer response 
from a multisine aileron input. 

A wingtip accelerometer response to a multisine aileron 
input is shown in figure 24. Figure 25 (top plot) shows 
the nature of a multisine in the wingtip response for a 35- 
to-5Hz sweep. Note the 13-15Hz pulse at the end of the 
main sweep (21 sec), characteristic of the multisine. The 


A TSVD-filtered reconstruction of the original wingtip 
response is in figure 27, which corresponds to the contour 
plot of figure 26. Both the SVD (bottom contour plot, 
figure 25) and TSVD contour and time history plots show 
dominance from the 12-13Hz mode just around 15 sec. 
Since the control system is active, this mode corresponds 
to the stabilator mode which is excited by the aileron 
multisine in closed loop. 
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Reconstructed Wingtip Response 
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Figure 27: Reconstructed TSVD-filtered F18-SRA wingtip 
accelerometer response from multisine aileron input. 


Time-Frequency Energy Density Distributions 
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Figure 28: Contour plots from original wavelet decom- 
position (top) and TSVD-filtered (bottom, log-scale), 
scalograms of F18-SRA fuselage lateral accelerometer. 


In a final example, figure 28 shows a persistent oscilla- 
tion of the wing bending mode at 8-9Hz (first fuselage 
and wing bending modes) from the fuselage lateral ac- 
celerometer for the same multisine aileron input. The 
top plot is the original energy density, and the bottom is 
the log-scaled TSVD-filtered distribution. The SVD dis- 
tribution is similar to that in figure 25, so is not shown. 
The TSVD defines the persistent mode very well while 
also detecting the dominant part of the chirp near 16Hz 
(10-15 sec). The final part of the energy density dis- 
tribution after 20 sec at 20Hz, however, is not detected 
since it is relatively less significant than the 8-9Hz modes. 
Time-dependent modal frequency and damping ratios es- 
timated with the TSVD filter are displayed in figure 29 to 
show the relatively low damping of the steady oscillation 
at 8-9Hz compared to modes at 19, 17, and 15Hz (10-15 
sec) corresponding to fuselage second bending, wing sec 
ond bending, and wing first torsion modes, respectively. 



Time, sec Timfi 


Figure 29: TSVD estimated modal frequencies (left) and 
damping ratios (right) from F18-SRA fuselage accelerometer. 


Conclusions 

Nonstationary dynamics data analysis is discussed from 
the perspective of time-frequency and multiscale resolu- 
tion analyses. Wavelet decomposition of dynamics data 
using a discrete harmonic basis constructed from the 
Morlet wavelet establishes a time-frequency representa- 
tion for modal parameter estimation and system identifi- 
cation procedures. Automated procedures to discern and 
extract regions of the time-frequency energy density dis- 
tribution that are significant for dynamics analysis are 
presented. The methods are derived from the singular 
value decomposition (SVD) of a time-frequency distn- 
bution. 

Examples of simulated sinusoids, a testbed aeroelastic 
system, aeroelastic flight data from a flutter experi- 
ment, and aeroservoelastic data from an F18 are used to 
demonstrate and compare the two SVD filtering meth- 
ods. One method using a standard SVD demonstrates 
the ability to discriminate the most significant and other 
less dominant dynamics. Its performance may depend 
on a time-frequency window determined by the input. 
Regions outside of the input time-frequency path relate 
to information about noise, unmodeled dynamics, and 
nonlinearity. 

Another method called the transformed-SVD, or TSVD, 
optimally rotates the singular vectors to create a more 
precise filter for finer tracking, but may ignore the less 
significant dynamics for the same order filter as the stan- 
dard SVD. Both methods complement each other and 
show promise for automatic feature extraction supple- 
mented with noise removal, residual analysis, and uncer- 
tainty estimation. These filtering methods contribute to 
nonstationaxy dynamics data analysis at multiple scales 
in the time-frequency framework. 
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